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Abstract 

We examine the linear stability of fluid interfaces subjected to a shear flow. 
Our main object is to generalize previous work to arbitrary Atwood number, 
and to allow for surface tension and weak compressibility. The motivation 
derives from instances in astrophysical systems where mixing across mate- 
rial interfaces driven by shear flows may significantly affect the dynamical 

evolution of these systems. 
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I. INTRODUCTION 



The stability of fluid interfaces in the presence of shear flows has been studied for almost 
half a century; and was largely motivated by the problem of accounting for observations of 
surface water waves in the presence of winds. As early as the 1950's, it was realized that 
classical Kelvin-Helmholtz instability could not account for the observed water waves 
(cf. [^,0), and efforts were initiated to study the full range of possible unstable modes by 
which interfaces such as represented by the water/air interface could become unstable. By 
the early 1960 's, the basic mechanism was understood, largely on the basis of work by Miles 
|pT|^3| and Howard [l^: They discovered that interface waves for which gravity provided 



the restoring force (e.g., waves that can be identified with so-called deep water waves) can 
be driven unstable via a resonant interaction with the ambient wind; this work was also one 
of the first applications in which resonant (or critical) layers played an essential role in both 
the physics and the mathematics. Work carried out at that time showed that the precise 
form of the vertical wind shear profile was critical to the nature of the instability; typically, 
it was assumed that the wind immediately above the water surface could be characterized 
by a logarithmic profile of the form 

U{z) = Uo + UiHz/6 + l) , (1.1) 

where Uo is the velocity jump (if any) at the water/air interface, z is the vertical coordinate 
(with z = marking the initial water-air interface), and 6 is the characteristic scale length 
of the shear flow in the air.Q The idea was then to demonstrate that surface gravity waves 



whose phase speed is given by c = \^gXA {g the gravitational acceleration, A the perturbation 
mode wavelength, and A = {p2 — pi)/(p2 + Pi) the Atwood number for a density interface 
between fluids of density pi [upper fluid] and p2 [lower fluid], with pi < P2) can couple to this 



^Such velocity profiles are commonly observed in the boundary layer of winds blowing over the 
surface of extensive bodies of water; cf. Miles (1957). 
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wind profile at a height z where c ~ Vi^z). At the time, it was not possible to construct a 
self-consistent description of the problem, such that a logarithmic wind profile automatically 
emerged from the analysis; and much of subsequent work has focused on establishing the 
nature of this wind shear profile (e.g., |^). Finally, we note that these studies have since 
been applied to a number of other contexts, including especially shear flows in atmospheric 
boundary layers, where they have been extensively expanded, including into the weakly 
compressible regime 0. 

Our own work is originally motivated by an astrophysical problem in which mixing at a 
material interface between two fluids with different densities is essential to the evolution of 
the astrophysical problem. Specifically, consider a white dwarf star, whose composition is 
almost completely dominated by carbon and oxygen. If such a star is in a close binary orbit 
with a normal main sequence star, then it has been known for some time (e.g., [^]) that 



accretion of matter from the normal star (largely in the form of hydrogen and helium) can 
lead to a build-up of an accretion envelope on the white dwarf which is capable of initiating 
nuclear hydrogen "burning" . This burning process can lead to a nuclear runaway, in which 
the energy released as a result of these nuclear fusion reaction is sufficient to expel a large 
fraction of the accreted matter in the form of a shell; such a runway is referred to in the 
astronomical literature as a "nova". The key element relevant to our present discussion 
is then that observations show that approximately 30% by mass of the ejecta are in the 
form of C+0 nuclei: since neither carbon nor oxygen are products of hydrogen burning in 
the accreted envelope, it must be the case that some sort of mixing process brought large 
amounts of stellar (i.e., white dwarf) carbon and oxygen into the overlying accreted material 
before envelope ejection. Furthermore, detailed analysis of the energetics of the runaway 
process has shown that simple hydrogen burning in the envelope cannot provide enough 
energy to power the observed nova; thus, additional energy release via "catalytic" nuclear 
reactions in which C+0 play important roles is required in order to match the observations 
(cf. [P^ p5| , pil| ,p5[]). Thus, from the perspective of both observed abundances of nova ejecta 
and consideration of the nova energetics, efficient mixing at the star/envelope interface is 



called for. Several possible mixing processes have been discussed in the literature, including 
undershoot driven by thermal convection in the burning envelope and Kelvin-Helmholtz 
instability; but detailed studies shown all of them to be ineffective in producing the required 



mixing (e.g., ||T6|,|T5|,|T^ ) . In this regard, the current astrophysical situation resembles the 
problem encountered by oceanographers in the 1950 's, as they tried to explain the observed 
mixing between the sea water-atmosphere interface. A new instability is needed to account 
for the observed mixing |]20[| . 



Following the previous oceanographic work, we explore the possibility that a critical- 
layer instability related to the coupling of stellar surface gravity waves to a shear flow in 
the hydrogen envelope - can account for the enhanced mixing rate. Thus, in this paper we 
embark on a systematic study of such an instability and apply our results to the specific case 
of mixing of C/0 to H/He envelope of white dwarf stars [^]. We note that similar scenarios 
can arise in a variety of other astrophysical systems, such as in the boundary layer between 
an accretion disk and a compact star, where mixing between fluids of different densities - 
as in the nova problem - is expected to play an important role. However, the earlier non 
astrophysical work largely focused on the case of very large density differences between the 
two fluids separated by an interface, and primarily considered the fully incompressible case 
(the weakly compressible case has been considered by 0). For the astrophysical case, the 
density ratios can be of order unity (for the nova case a typical value would be Pi/p2 ~ V^O) 
and the Mach number for the interface between the accretion flow and the white dwarfs 
surface can range from very subsonic to of order 0.2. The aim of this paper therefore is to 
extend the shear flow analysis to arbitrary density ratios , shear and compressibility. We 
provide estimates of the growth rates of unstable surface waves, and determine the regions 
in the control parameter space that correspond to different instabilities for different physical 
situations. 

This paper is structured as follows: In the next section we define the problem to be solved 
more precisely; §fT^ describes the linear analysis for the incompressible, two-layer case. § P^ 
and §0 describe, respectively, the inclusion of surface tension and extension to compressible 



flow of low Mach numbers. We discuss and summarize our results in the final section. 

II. FORMULATION 

The starting point of our formulation is the identification of the appropriate material 
equations of motion. This issue has been well-discussed in the literature, including the mo- 
tivating white dwarf case 0: in general, we can expect the gaseous surface and atmosphere 
of such stars to be well described by the single fluid equations for an ideal gas. More specifi- 
cally, the length scales of the physical mixing processes discussed here are all far larger than 
the Debye length, so that the ionized stellar material can be considered to be neutral; and as 
long as the stellar magnetic fields are weak (e.g., (3 = gas preassure/magnetic preassure ^ 1) 
we can ignore the Lorentz force. Furthermore, the ratio of the spatial scales of interest to 
the Kolmogorov scale is large (typically > 10^) so that we are in the large Reynolds number 
limit, and viscous effects on the motions of interest will be negligible. As a consequence, the 
Euler equations will be describing our system. 

We consider a two-dimensional flow with x the horizontal direction and y the vertical. 
The system consists a layer of light fluid (density pi) on top of a layer of heavy fluid 
(density p2)- In most of our analysis pi and p2 are constant in each layer, and in the 
most general scenario both layers can be stratified (densities are functions of y). The two 
layers are separated by an interface given hj y = h{x;t), which initially is taken to be flat 
{y = h{x;0) = 0). The upper layer (pi) is moving with velocity U(?/) in the x direction 
parallel to the initial flat interface, while the lower layer (p2) remains still. 

As already mentioned, the instability of such stratifled shear flow has been investigated 
(cf. p0| , p!3| ), albeit under limited physical circumstances. We study this problem in full 
generality, allowing for a variety of effects (including broad ranges in the values of the 
Atwood number /gravity and in compressibility) with the motivation that one can establish 
the role of the relevant instabilities under more general astrophysical circumstances than the 
restricted case of nova-related mixing which motivated our study. 
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A. The general problem 



A wind (shear flow) is assumed to flow only in the layer of light fluid (pi) and is zero in 
the heavy fluid (^2)- Within each layer, the governing equations are the continuity equation 

dtp + W ■{pu) = Q, (2.1) 

and the two-dimensional Euler equation 

pdtu + pu-Vu = -VP + pg. (2.2) 

The equation of state closes the system, which is expressed in dynamical terms: 

where 7 is the polytropic exponent. The background density and pressure are in hydrostatic 
equilibrium, dyPo = —PoQ- The basic state is then defined by a shear flow (U(?/)) in the 
upper layer, and hydrostatic pressure (Pq) and density profiles (po)- We perturb around this 
basic state 

u = \J{y)x + u', p = p^{y)+p\ P = Po(y) + p', (2.4) 



and study the growth of the perturbations (primed variables). From equation p.3| , the density 
and pressure perturbations satisfy the relation 

= Cs-^ + w [gpo + c,pay), (2.5) 

where w' is the vertical component of the perturbation velocity, g is the gravitational accel- 



eration, and Cs = y'jPo/Po is the sound speed for the background state. Upon expanding 
the perturbations in normal modes e^'^i^-'^^)^ obtain the linearized equations in the per- 
turbation quantities (where we have dropped the primes for convenience): 

ik(\J — c)u +wdyU = —ikp~^p, 

ik{U - c)w = -p-^dyp - p-^gp , 

(2.6) 

ifc(U — c)p = —poiiku + dyw) — wdyPo , 

z/c(U — c)p —wgpo = cj. {ik{U — c)p + wdyPo) . 
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The above equations form an eigenvalue problem for the complex number c. One imme- 
diately sees that the incompressibility condition V ■ m = can be obtained by taking the limit 
Cs — oo. Our problem simplifies greatly with this assumption. Therefore we first present 
our results for the incompressible case, and then examine how compressibility modifies the 
stability properties. 



B. The incompressible case 

For the incompressible case we define a stream function such that u = dy"^ and 
w = —dx"^- The 2-D Euler equation thus reads 

SiV^^ - ^.V^^, + ^,V'^, = 0, y^h. (2.7) 

The total stream function \E' = + ip consists of a background stream function \l/o = 
/o \J{z)dz and a perturbation ip = (f){y)e'^''^^~'^^\ The linear equation for is the well-studied 
Rayleigh equation, 

0"-(fc'+U^)0 = O. (2.8) 

The boundary conditions at the interface for the continuity of the normal component of the 
velocity and pressure are: 

(U - c)h - (f)^ = 0, (2.9) 
A[p,((U - c)0' - U»] + gh{p, - p2) = 0, (2.10) 

where A indicates the difference across the interface, and h is the amplitude of the perturbed 
interface, h = he^H^-(*)_ 



C. The compressible case 

For the compressible case, where Cg is comparable to the background shear flow and the 
density stratification is non-negligible on the scales of interests, we start from the full set of 



equations |2]^. We obtain the following equations by eliminating p and u: 



(2.11) 



Po{k'^U^ + gkg + gks)w = ikUoidy + kg)p , 
tk\Ullcl - l)p = PokikgUa + dyUa - UGdy)w , 

where Uq = f/ — c is the Gahlean-transformed velocity in the reference frame of the wave, 
ks = Po^dyPo is the inverse stratification length scale, and kg = g/c1. We further simplify 
the equations by applying the transformation 

p = f~^p, w = iUcqf, po = Pof with / = e-^o ^a<y^y\ 



Equations |2.11| are then rewritten in terms of these new variables as follows: 

p{k'^Ul + gkg + gks)q = kdyP, 
k\l-Ul/cl)p =PokUldyq, 

which can be combined to give 

poUcdyq 



(2.12) 



Po 



k'U^ + g{h + kg) q = 0. 



(2.13) 



We re-write the above equation into canonical form. The resulting equation is similar to the 
Rayleigh equation for the incompressible flow, except for an additional stratification term 

r^2 



-g{ks + kg)/\](.(l): 



_ kg -\~ k, 



9 _|_ "y 







(2.14) 



where = fc^ (1 — Uq/cI), Ug = kUc^fpo/i^ and (p = qVc/k = —ik ^y/fkiw (1 — Uq/c^) \ 
It can be shown that the stratified Rayleigh equation can be recovered by taking the limit 
of c. 



oo. 



k' + 9T^ + 

G 



ks ^'^{\]G^/p^) 







(2.15) 



Furthermore, we recover the unstratified Rayleigh equation in the same limit, if fc^ + = 
(which corresponds to an adiabatic atmosphere, as we will show later on). Finally, the 
boundary conditions at the interface are expressed in terms of Ug and (p: 



q 



US U5 



h 



(2.16) 



using the continuity of q and integrating |2.13| across the interface we obtain 



UGdy(t>-{dyiJG)<p-g~p<P/VG =0. (2.17) 



D. Wind profiles 

In general, it is not trivial to determine the wind profile: strictly speaking, the wind 
profile should be determined as part of the solution of the evolution equation for the wind 
shear interface. However, it has been customary to simplify the problem by assuming an 
a priori analytical form for the wind profile, and to use this in order to study the stability 
properties of the interface; thus. Miles |^ used a logarithmic wind profile from turbulent 
boundary layer theory to model the wind profile in the air over the ocean. In this example, 
the turbulence level in the wind is simply defined by the scale height of the wind profile, 
which in turn simply depends on how "rough" the boundary is. 

In our formulation we shall also assume the wind profile to be of simple form, and 
scale distance with respect to the length scale of the wind boundary layer. In order to 
explore the sensitivity of our results to the nature of this wind boundary layer, we will 
examine two different kinds of wind profiles: the first is the logarithmic wind profile U(?/) = 
Uq + Ui ln(a?/ + 1), which is derived from turbulent boundary layer theory for the average 
flow above the sea surface; the second is given by U(?/) = ?7i tanh(a?/), which has the more 
realistic feature of reaching a constant finite flow speed above the interface. 

III. LINEAR ANALYSIS: INCOMPRESSIBLE CASE 

We start with the stability analysis of the incompressible case with constant densities in 



the two layers. The fluid is described by the Rayleigh equation |2.8| within each layer; and 
we ignore surface tension for the time being. We solve the following equation in each layer: 



^yy 



{^' + ir^M = ^^ 0U±oo = o, (3.1) 



with boundary condition (at y = 0) 



p2kc^ - pi {U - cf(t)y -{U - c)Uy - g{p2 - pi) = 



(3.2) 



where we have normahzed (f) by setting <^|y=o 



1. 



We scale lengths by a 



, the characteristic length of the wind profile and the velocity 



by the reference velocity Ui. The dimensionless equation thus reads 



K - {k' + 



V -c 



)0 = O, 0|,=o = l, 01 



y=oo 0, 



(3.3) 



and the boundary condition at the interface now becomes 



KC^ -r {Vo- Cf(Py - {Vo - C)Vy - G{1 - r) = 0. 



(3.4) 



where C = c/Ui, K = k/a, G = g/Ula, Vq = U(0)/Ui, and r = pi/p2- 

For a given wind profile, the system then is characterized by the four parameters 
{K, G, Vo, r). Parameter G measures the ratio of potential energy associated with the surface 
wave to the kinetic energy in the wind. The Richardson number defined in stratified shear 
flow is not useful in quantifying the stability in our case. However, as will be shown later, 
we find parameter G to be a good substitute in describing the effect of stratification on the 
surface wave instability. In the case of accretion flow on the surface of a white dwarf G ~ 1, 
while in the case of oceanic waves driven by winds, 0.1 < G < 1.0. Table | lists the values 
of G for a variety of physical conditions. 

The aim of our linear analysis then is to find the value of G in the complex plane as 
a function of these 4 parameters, and to establish the stability boundaries in the space 
{K, G, r, Vo); note that in our convention, S^{G} > implies instabihty (where refers to 
taking the imaginary part). 



^In oceanography, such a length scale is referred to as the "roughness" of the wind profile. 
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A. Kelvin-Helmholtz modes and Critical Layer modes 



Before solving this problem, some general remarks about the set of equations O - O are 
required. We observe that in the inviscid limit, if C is an eigenvalue, then so is C*; therefore 
we will have a stable wave only if ^^{C} = 0. If that is the case, then at the height where 
Vc = V{ycr) = C (assuming such a height exists) the Rayleigh equation has a singularity; 
this location y = ycr is called the critical layer, and is well-discussed in the literature P,PD[] 

The existence of such a critical layer is crucial for the presense of instability. One can 
prove (Appendix A) that our system can be unstable only if Cr = 3?{C} < Vmax- For the 
case that Vo = there always exists a point in the flow where Cr = V < Vmax for all unstable 
modes. We denote this point as a critical layer even if C is complex i.e., Ci = '^{C} ^ 0; and 
thus there is no singularity. However, if Vo ^ 0, such a point might not exist (e.g. if C,. < V^). 
In that case the only mechanism that can destabilize the flow would be a Kelvin-Helmholtz 
instability. These two kinds of instabilities exhibit very different properties, both in terms of 
the physical mechanisms involved as well as in the mathematical treatment required. Hence 
we need to distinguish between (i) modes becoming unstable due to the discontinuity of the 
wind proflle (from now on called KH-modes), and (ii) modes becoming unstable due the 
presense of a critical layer (from now on called CL-modes). 

The KH-modes have been studied for over a century; here we summarize the results for 
a step function wind proflle and some of the features can also found for other wind proflles 
1^. The dispersion relation is given by 

P2c' + Pi(c + Uf = (P2 - Pi)9/k, (3.5) 

where U is the jump in velocity across the interface between pi and p2- The pressure 
perturbation pi{c + U)"^, providing the driving force for the instability, is always in phase 
with the wave, and is independent of the wavelength. The restoring force (p2 — Pi)g/k 
on the other hand is proportional to the wavelength, and so we have instability when the 
wavelength is sufficiently small for the pressure to overcome the restoring force. In more 
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physical terms, the flow stream lines above the crests of the perturbed interface wave are 
compressed, and above the troughs are decompressed. According to Bernoulli's equation 
the pressure above the crests is therefore decreased, and is increased above the troughs. The 
wave thus becomes unstable when these destabilizing pressure forces exceed the stabilizing 
effects of gravity. The dispersion relation 



-U± 



- A I - T^^.U^ (3.6) 

P1+P2 \^k\pi+p2j {pl+p2)^ 

also shows that the growth rate becomes positive only for wavenumbers k > g{p\ — 

pl)/iU'piP2). 

The CL-modes behave very differently. The solutions of near the critical layer for 



small or zero Ci have a singular behavior. The two Frobenius solutions at the point where 
y = Ucr are given by 

where z = y — ycr (subscript cr means "evaluated at the critical point"). The singular 
behavior appears in the first derivative of (f)},. The singularity is removed either because 
Cj 7^ 0, in which case the Frobenius solutions have the same form but ycr is now complex 
(so z never becomes zero); or because viscosity becomes important in this narrow region, 
in which case the inner solution can be expressed in terms of generalized Airy functions 
0. In either case, the basic result is that there is a phase change across the critical layer, 
by which we mean that if = a(f)a + ^'Ph is the solution for the stream function above 
the critical layer, then the solution below would be = (a + i7Tb)(f)a + in the previous 
formula. Physically this means that the perturbation wave above the critical layer is not in 
phase with the wave below this layer. Moreover, when we apply the boundary conditions at 
the interface, since dy(f)\o is now in general complex, the solution of equation will give a 
complex value of C. That is, the pressure gradient reaches minimum value not on top of the 
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crests, but rather in front of the crests, where gravity does not act as effectively as a restoring 
force. In particular, the destabilizing force is now non-zero at the nodes of the boundary 
displacement field (i.e., where h = 0), where the gravitational restoring force vanishes, but 
where the vertical velocity of the interface is maximum; thus, the forcing resembles pushing 
a pendulum at its point of maximum velocity but minimum displacement. Note that in this 
case, there is thus no cut-off for CL-modes corresponding to the wavenumber cut-off due to 
gravity for KH-modes. 

Having discussed the physical mechanisms for destabilization, we now turn to the im- 
plications for our choices of initial wind profiles. For wind profiles with Vo = one notices 
that if we assume (py to be real and known, then the complex eigenvalue C is obtained by 



solving equation 3.4 



r ± + 4G(1 - r){K ~ r^y) 
^ ^ 2{K - r<f)y) ' 

which will have a nonzero imaginary component only if (py is positive and {K — rcpy) < 0. 
However, the negative real part of C implies that the surface wave would be traveling in 
the direction opposite to that of the wind — this case can be excluded on physical grounds 
(a more rigorous proof is given in appendix A where we show that Cr > 0). Thus, the 
mechanism that gives rise to the unstable KH-modes can be excluded. Thus we conclude 
that surface waves become unstable in this case only if a critical layer exists. If we use the 
logarithmic wind profile, we obtain unstable waves for all wavenumbers because ln{y + 1) is 
an unbounded function, therefore a point y where Cr = V{y) exists for every value of Cr- 
This however is not true for the tanh wind profile. Because waves with C,. > Vmax are stable 
and Cr (in the absence of surface tension) is a decreasing function of K, there must be a 
lower bound on K, Kmm, so that waves with K < Kmin are stable, and unstable otherwise. 
The value of Kmin in general will depend on the exact form of the wind profile. In appendix 
B we find the exact value of Kmin for a wind profile of the form V = 1 — e^^. 



G{l-r)+r- r^/(G(l - r) + r)2 + (1 - 
T 



Kmin 1 -„2 ' ('^■'^) 
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We remark the following about the previous formula. First of all we note that although 
the previous result holds only for the specific wind profile used, it can provide a general 
estimate of Kmin- Moreover we note that, unlike the Kelvin-Helmholtz case, in the limit 
r — > 0, Kmin remains finite and equal to G (however, the growth rate goes to linearly with 
r i.e., KCi ~ r); this confirms that for small density ratios CL-modes dominate. Finally 
by writing the wind profile in its dimensional form U = Ui{l — e~"^) and taking the limit 
a ^ oo (which takes the wind profile to the limiting form of a step-function, U = Ui for 
y > and Ui = otherwise) we get kmin = — r)/U'^ which is different from the result 
Kelvin-Helmholtz instability gives. We therefore conclude that different limiting procedures 
lead to different results. 

The situation is more complicated if > 0. No critical layer exists for Cr smaller than 
Vo- Hence, modes of sufficiently large K become stable (e.g., there is a upper bound on K for 



the unstable CL-modes). One notes that C = Vo = JG{1 — r)/i^ is a solution to equation 



3^ . This solution corresponds to the case where the critical layer is right at the interface. 
For slower modes than this (C^ < Vo) a critical layer will not exist, and therefore the surface 
gravity modes will be stable. Thus CL-unstable modes exist only for K < G{1 — t)/V^^; 
this result has been confirmed numerically. As K is increased, the discontinuity of the wind 
profile becomes important and Kelvin-Helmholtz instability rises. The system therefore will 
be unstable for K < Kcl and for K > Kkh, where Kcl is the upper bound of the CL- 
modes and Kkh is the lower bound for the KH-modes. This implies that there is a band 
of wave numbers Kql < K < Kkh that corresponds to stable modes, and separates the 
two unstable wavenumber domains. However, as we will show later, for some values of the 
control parameters this stable region disappears, and the two instabihties overlap. 

B. Small density ratio 

We are now ready to present results from the linear analysis for the logarithmic and the 
tanh wind profiles. The existing literature has primarily covered the case of small r, with 
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the other parameters assumed to be of order one. In contrast, we are interested in covering 
a wider range of the control parameters, and thus provide a complete description of the full 
dispersion relation C = C{K). We therefore briefly summarize Miles' results and move on 
to the general case. 

Assuming the mass density ratio r is a small number (which is true for the air over water 
case) and the other parameters are of order one. Miles [^] expanded the eigenfunction and 
the wave velocity C with respect to r 

(j) = r(j)o + r^(j)i + r^(j)2 + h.o.t., C = Cq + rCi + r^Ci + h.o.t. ; (3.10) 

one then obtains the zeroth order solution as a linear gravity wave with constant amplitude 



and phase speed Cq = ^G/K. At first order 0(r), one finds 

2A'C„Ci - (l'„ - C„)%,l,„ + (V„ - C„) + G = 0. (3.12) 

The asymptotic expansion breaks down at the critical point y = Ucr since to first order Cq 

is real. Using the phase change of iir rule across the critical layer from theory [Q, Miles 
obtains the growth rate of the perturbation at leading order in r: 

^^^'^ = 2^^^° ~ go)Mg.0o} = {%] \<Pcr\\ (3.13) 

\ y / cr 



where the last relation is obtained by multiplying equation p.3| with the complex conjugate 
of and taking the principal value integral, with the contour going below the singularity; 
the subscript "cr" means evaluated at the critical point. 

The first case we examine is when the velocity at the interface is zero. This simpli- 
fies things slightly because, as we discussed before, there are no Kelvin-Helmholtz unstable 
modes in this case. The dispersion relation '^{C{K)) is shown in figure |l|a,b for the logarith- 
mic and for the tanh wind profile for various values of G. The only difference between the 
two wind profiles appears at small wavenumbers: the tanh wind profile (whose asymptotic 
wind speed is bounded) does not permit waves traveling faster than the wind to become 
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unstable. For this reason there is a cut-off which can be found in our small r approximation 
to be at = G for the tanh wind profile. 

Next we look at the case where Vq > 0. Now we have both modes present. As discussed 
before, the CL-modes are stable for wavenumber K < Kcl = ~^)/^o- The KH-modes 
will appear when VoK/G increases to order 1/r. If we denote by K^h the minimum value of 
K that the KH instability is allowed, then for the KH modes in the small r limit K scales as 
1/r, and one can perform a regular perturbation expansion for large K, small r (Appendix 
C) to find: 



1 d'^V 1 r cfiV d'^V d V 



K-^ + ..., (3.14) 



and 



KKH = ^--iVoV'\o + G) + ... . 

Vo^ r 

The above resembles result for a step function wind profile except for small corrections 
due to the non-constant velocity profile. Thus for small density ratio, the difference between 
the two modes is as discussed in § P^1I A| . We will discuss the two instabilities in more detail 
in the r = 0(1) case. 

C. Large density ratio 

For large density ratio, we solve the system of equations IHT^-P]^ directly. We focus on the 



instability properties of special interest, such as the maximum growth rate, the wavelength 
of the fastest growing mode, and the dependence of the the stability boundaries on the 
parameters of the model. First we present results for cases where the wind has zero velocity 
at the interface (Vo = 0) in figures ||a-c and |^-c. We solve equations ^73| - ^]4| numerically 



using a Newton- Ralphson method; the results for both wind profiles logarithmic and tanh are 
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presented together for comparison. The plots sugest that for small enough r the dependence 
on r is linear (e.g. the r = 0.001 case is proportional to the r = 0.01 case by exactly a factor 
of 10.0). For larger values of r the dependence is stronger than linear, and the smaller K 
modes seem to become more unstable. 

We have repeated these calculations for the case Vo ^ 0; the results for the 55(C) are 
shown in figures ^-c and |^a-d. In this case we have to distinguish again the two different 
kinds of instabilities. The distinguishing factor for the Kelvin-Helmholtz instability (most 
prominent in the discontinuous wind profile) is that the growth rate is positive only for 
wavenumbers larger than some lower bound. However, the critical-layer instability, which 
owes its presence to the phase change in the critical layer, has an upper bound in wavenumber 
for instability. Thus in general there exists a band of wavenumbers for which both modes 
are stable. The difference between small and large density ratios is that the two instabilities 
are not necessarily separate, and they do overlap for some parameters. The criterion for a 
critical layer to exist in this case, \JG{1 — r)/ K > Vo, provides a upper bound on K for 
unstable CL modes. An exact solution for the upper boundary is not available, but the 
asymptotic behavior of the second boundary, for large K and for small Vq (cf., equation 
3.141 ) suggests that it takes the form of K~^^'^; thus, the two stability boundaries are not 
expected to cross in the large wavenumber limit. However, the two boundaries do meet for 
small K and large Vo, as it can be seen in the stability boundary plots ^-c. 

D. General features of the CL instability 

The main goal of this paper is to establish the relevance of the critical-layer instability 
under various astrophysical or geophysical conditions. Results from our linear analysis allow 
us to identify the most unstable modes in different parameter regimes (and thus physical 
situations). Furthermore, the maximum growth rates give us an estimated time scale of the 
nonlinear evolution, and the length scale of the fastest growing mode allow us to estimate the 
thickness of the mixing layer as instability grows; this allows us to provide rough predictions 
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of the physical conditions for which more exact fully-nonhnear calculations should be carried 
out to establish the mixing zone properties more precisely. With this motivation in mind, 
we show here how these two quantities behave as functions of the physical parameters G, r 
and a. 

In figures l^a-b we have plotted the maximum growth rate 7 = '^{CK}max of the pertur- 
bation as a function of the control parameter G for the two wind profiles and for different 
values of r. It is clear in all cases that there is an exponential dependence on G for G ^ 1. 
This might be expected because increasing gravity leads to an increase of the real part of 
C; therefore the imaginary part of (p, that falls exponentially with the distance from the 
critical layer, will have an exponentially smaller component at the interface. Furthermore, 
as the density ratio r approaches unity, the dependence on gravity becomes weaker. We 
plot the slopes of the curves from figures ^ as a function of r in figure The dependence 
on r is roughly linear (deviations from linearity will not be important since r only takes 
values in the range < r < 1). This allows us to write an empirical scaling law law for the 
dependence of the growth rate on the control parameters: 

7^„, = ^ /?re-"(i-'^)^ (3.15) 

For the logarithmic wind profile we found a ~ 2.8 and (3 ~ 0.10; while for the tanh wind 
profile we found a ~ 2.9 and j3 ~ 0.18. 

In figures p!0|a,b we have plotted the wavenumber of the most unstable mode (whose 
growth rate corresponds to equation |3.15|) as a function of G. We see that the wavelength 
of the most unstable mode highly depends on the wind profile length scale. In particular, 
for fixed density ratio r, Xmax ~ the dependence on gravity or on r is weaker. 
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IV. SURFACE TENSION 



For the sake of completeness, we have also examined the case in which surface tension at 
the density interface is included^. We again assume a wind shear profile of the form ln(?/ + 1) 
and tanh(y). The only change in our set of equations to solve is then an additional term in 



the boundary condition, equation |3]J. Hence 



KC^ -r 



{Vo - CY<Py - (K, - C)Vy - G(l - r) - TK 



0, 



(4.1) 



where T = aa/ {p2V\)) and a is the surface tension (a = B"^ /(2njiK) for the case of the 
magnetic field [§]). We show the resulting solutions, namely the dispersion relations, in fig- 



ures |TI|a-c and |l2|a-c. As expected tension decreases the growth rate and becomes important 
only in large wavenumbers. 

An important result, which we have not previously seen derived, is that in the small 
density ratio limit, the real part of C (to zeroth order in r) is Cq = JG/K + TK which has 



its minimum value C^m = v^(G'T)^/'* at K = y G/T. Thus, for the case of a bounded wind 
profile (such as the tanh profile) there is a minimum value of Ui, given by Gmin, so that a 
critical layer can exist. We remind the reader that a similar minimum velocity bound also 
exists for the Kelvin-Helmholtz instability, and is given by 



U > 



\ 



2 



650 cm/sec , 



where we have retained only terms of first order in r. For the CL case we have instead 



U > 



\ 



2J — ^ 20 cm/sec 
V P2 



which differs from the previous bound by a factor of ^/r. (The numerical values shown here 
are derived for the case of air blowing over water.) This illustrates the fact that for low wind 
conditions, the CL instability dominates the KH instability for driving water surface waves. 



^ We note that a magnetic field in the lower fluid, whose direction is aligned with the flow, would 
lead to an equivalent treatment (see for example Q §106). 
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V. COMPRESSIBLE CASE 



Finally we consider the compressible case. We will consider a compressible fluid in the 
upper layer with sound speed Cs{y) and an incompressible fluid below. The dimensionless 
equations we have to solve now are: 



p A', + A'^ 



P\y=0+ 



T/2 



Vg 







(5.1) 



V^\y=0dy^-VG\y=0dyVG\y=0 - G{l - t) ^ ^ 



(5.2) 



where 



and 



dyP 



ap 



KVG^P/P\y=^+ 

Vg = — with r = p\y=Q+ / p\y=Q- . 



K 



We will assume for simplicity an adiabatic atmosphere, 



p = p\y=o+ f 1 - (7 - ^)^9y 



1 

1-7 



(5.3) 



(5.4) 



This assumption, which is commonly used in the atmospherical sciences to simplify the 
physics involved, has the advantage that Kg + Kg = 0, so our equation becomes by one 
order less singular, and therefore becomes easier to solve. 

We will not deal here with supersonic flows since in most astrophysical realms in which 
interfacial wave generation plays an important role (viz., on white dwarf surfaces) the rele- 
vant flows are thought to be subsonic; for this reason we will consider only the tanh wind 
proflle. We shall also deal with small values of G, so that the pressure scale height is large 
and the breakdown of the adiabatic assumption at values of |/ ~ will not affect us 
either. 
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The dispersion relation for different values of Cs and for a tanh wind profile is given in 
figures |13|a,b. Compressibility, as it can be seen from the figures, decreases the growth rate. 
This is an expected result, since our system has now more degrees of freedom (e.g., now 
the perturbation stores thermal energy as well). We conclude, however, that the deviation 
from the incompressible case is not very large, even for relatively strong (but still subsonic) 
winds. 



VI. DISCUSSION AND CONCLUSION 

In this paper we have explored the linear instability properties of wind shear layers in the 
presence of gravitational stratification. Our principal aim was to explore the full parameter 
space of the solutions, defined by the four parameters K (the perturbation wavenumber), 
G (related to the Richardson number, and measuring the relative energy contributions of 
the gravitational stratification and the wind), Vo (the velocity discontinuity at the density 
interface), and r (the density ratio). 

We have distinguished between the two different kind of modes (Kelvin-Helmholtz modes 
and Critical Layer modes) existing in our model and constructed stability boundaries for 
those, as well as the dependence of these boundaries on the given parameters. As we discuss 
later, the non-linear development of the instability (and therefore mixing) will crucially 
depend on the kind of modes that become unstable; therefore an investigation of the stability 
boundaries is necessary before the study of the nonlinear regime. An important result 
also derived from our analysis, allowing us to make predictions on the importance of the 
instability and on the nonlinear development, is the scaling of the growth rate with the 



parameters G and r in § [111 D| ; our results show that for (1 — r)G ^ 1 strong mixing is not 
expected. This result has an interesting implication: although, as expected, strong gravity 
(e.g. G) inhibits mixing, one might still experience strong instability in the large G case if 
the density ratio of the interface is close to but not equal to unity. Finally we investigated 
the effects of surface tension and compressibility. With the inclusion of surface tension, we 
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obtained a lower bound on Umax for the instability to exist. We also found that for subsonic 
winds the instability growth rate weakly depends on the Mach number. 

As we have shown, there are significant differences between the CL and the KH modes, 
both in the parameter ranges in which the instability can occur (e.g., the stability bound- 
aries) and in the magnitude of the growth rate; these differences can be expected to result 
in different nonlinear evolution of the underlying physical system. For example, it is well- 
known that CL instability in the air over water case is responsible for generating waves for 
winds of magnitude below the threshold for which Kelvin-Helmholtz instability exists ^0 



An important aspect not discussed as yet is the case in which the spatial density vari- 
ation is smooth instead of discontinuous. In our simplified model of a sharp interface, the 
distinction between CL modes and KH modes emerges naturally from our analysis, simply 
based on the existence or absence of a critical layer. In the more realistic (astro)physical 
case, however, sharp velocity and density gradients do not exist. For this reason, we need 
to generalize our definitions for the two kinds of modes. We proceed by considering the 
physical mechanisms involved in the instabilities: In the KH mentioned above, the 

pressure perturbations are in phase with the gravity wave amplitude; and the wave becomes 
unstable when pressure overcomes the restoring force (e.g., gravity). An immediate con- 
sequence of this is that when the restoring force is overcome, it no longer plays a role in 
the wave propagation, and therefore the real part of c is independent of the restoring force, 
i.e., independent of g. This argument can be confirmed by examining the results for the 
step-function wind profile, where we see precisely the predicted behavior. 

In the CL case, it is instead the out-of-phase pressure component that drives the wave 
unstable; in this case, the destabilizing pressure force does not strongly modify the restoring 
force (here again gravity). Hence the real part of c is only weakly modified when the mode 
becomes unstable, and therefore the wave continues to propagate with its "natural" speed 
while going unstable. These properties of wave destabilization, which affect the dependence 
of the real part of c on the restoring force, can therefore be used to distinguish between the 
CL and KH modes. Thus, in the more general case, we shall refer to the modes that become 
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unstable due to an in-phase pressure component as KH modes; their propagation speed is 
independent of (or at most weakly dependent on) the restoring force. In contrast, we shall 
refer to the modes that become unstable due to an out-of-phase pressure component as the 
"resonant modes" (since the name "critical layer" is not appropriate for the general case); 
their propagation speed depends strongly on the restoring force.Q 

By extrapolating our results to the smooth density variation case we conclude that KH 
modes are likely to appear in cases in which the inflection point of the wind profile, or the 
region in space in which U changes, is at the same height as the region which the density 
changes. We note that the KH instability, as defined in § [111 A| , is a limiting case of such 



wind profiles. In contrast, resonant modes are more likely to appear when the regions of 
velocity and density change are well separated, where the coupling between an existing 
gravity wave and a critical layer above can lead to a 'resonant' behavior as described above. 
This expectation is supported even further by the observation that in the case of a smoothly 
stratified fluid the stratification term becomes dominant in the critical layer and the phase 
change is no longer in but depends on the Richardson number. 

We also note that, in the past literature on shear flow instability, much attention is 
focused on the KH instability according to our definition. For example, models with U ~ 
tanh(?/) and dy ln(p) ~ tanh(?/) or p ~ exp(— ?/) have been studied in the linear, weakly 
nonlinear, and fully nonlinear regimes P,|I|,|],|l2| , but they all fall into the KH category (as 
we defined it); complete study of resonant modes, though, has not been fully investigated. 
This is a gap we seek to close in our further work. 

From our results and the ranges of our physical parameters in Table H we can estimate 
the growth rates of the wind driven waves as well as their typical wave length. For the 
astrophysical problem we are interested in, we conclude that the growth rate can be as large 



^The words 'weakly' and 'strongly' are used here because it is expected that there is a smooth 
transition from the one case to the other as we change our physical parameters. 
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as 10~^s~^ with typical wave lengths of the order of a~^; these results were obtained using 
the empirical formulae p.l5| and p.9| . For our motivating astrophysical application, i.e., the 
nova mixing problem, the results shown in |^ are especially important. First, we note that 
the interface between the stellar surface (at the typical density p2 = 3800gcm~^) and the 
accreted envelope (at the typical density pi = AOOgcm'^) is a material (gaseous) boundary 
at which one would not expect any free slip. Thus, we would expect Vo = Uo/Ui (in, for 
example |1 . 1| ) to be very small, and essentially zero. Consider then Panel(b) of ^ (for which r 
takes on the astrophysical relevant value): we see that for small V^, the interface instability 
is completely dominated by the resonantly-driven modes; the classical KH instability only 
appears at very large wave numbers, e.g. very short wavelengths and therefore is unlikely 
to matter in the nova mixing problem. To go beyond this will require further investigation 
of the nonlinear evolution of the CL instability and is currently under investigation; more 
information on the astrophysical model, is provided in |2^. We discuss some of the issues 
relevant to the nonlinear regime next. 

Finally, we comment on the possible nonlinear development of the two types of modes we 
have studied. The nonlinear evolution of a KH-mode is in fact well-studied in the literature 
|I|J^|T^ and is known to lead to a mixed region of width roughly equal to the wavelength of the 
mode; indeed, mixing proceeds in this case until (roughly speaking) the Richardson criterion 
holds throughout the flow. In the case of CL modes, the nonlinear evolution is affected by 
the fact that a new length scale enters the problem, namely the height of the critical layer, 
Uci which can be substantially larger than the mode wavelength. Thus, one might expect 
that mixing proceeds until heights of order are reached by the mixing layer, and therefore 
we expect more extensive mixing. Clearly, the next steps in this study involve investigation 
of the weakly nonlinear regime (to examine supercriticality and possible saturation of the 
modes), as well as the fully nonlinear regime (through numerical simulations). A particularly 
interesting question is to what extent the expected wave breaking of the CL modes (cf. p7[ ) 
can lead to enhanced mixing at the shear /density interface. 
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APPENDIX A: EXTENSION OF HOWARD'S SEMICIRCLE THEOREM 

We begin from 




J/=oo 



= 



(Al) 



and 



(A2) 



where G is the restoring force (G = G(l — r) for the simplest case), < G, < K, and we 
assume Cj 7^ 0. Let Vq — V — C, and let ip — (I)/Vg and D = dy. Note that 



V^Dcf) - VG(f)DVG = V^DjP . 



The boundary condition can then be written as 



rV^D^Ij = KC^ - G . 



(A3) 



Prom (Al) we obtain 



(A4) 



multiplying the last relation with Voij^* and integrating, we obtain: 



so that 
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D[i;V^D,l;] - V^\D,I;\' - V^K'\i^\' = 



and 



using the normalization condition, and denoting by Q{y) — \\D'^\ + X^l-^p] > 0, we then 
have 

1 ..oKC^-G 



SO that 



KC^-G 



--rl {V-CfQdy. 



Taking the imaginary part, we obtain 

K{2CrCi) 



roo 

■ 2Ci{V-Cr)Qdy 
Jo 



(A5) 



KCr 



{V - Cr)Qdy 



Therefore, 



(A6) 



0<Cr<Vr 



r ^ " max 



(A7) 



i.e., a wind cannot generate waves travehng faster than its maximum speed. Now, taking 
the real part, we obtain 



K{C^ -Cf)-G 



-r / [V - 2VCr + - Cf]Qdy , 
Jo 



(A8) 



or 



K(C^ - Cf) - G 



roo roo roo 

/ V^Qdy - 2Cr / VQdy + (C," - Cf) / Qdy 
Jo Jo Jo 



or 
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ir(C2 - ) - G 



V'Qdy - 2Cr I + Cr Qdy\ + {C^ - C}) Qdy 



so that 



< 



K\C\^ + G 



r [V'-\C\']Qdy 
Jo 



and 



(A9) 



(AlO) 



which is the sought-for result. 



APPENDIX B: LOWER BOUND ON THE CL-UNSTABLE MODES 



Consider the wind profile V = 1 — e ^; then the equations 3.3 and 3.4 become 



"yy 



1 - e-^ - C 



(Bl) 



and 



- r[C^<py + C] - ^(1 - r) = 0. 



(B2) 



We are interested in the value of K for which our system becomes marginaly unstable. From 
the extension of Howard's semicircle theorem to our case we know that for Cr greater than 
Vmax the system is stable, so the instability is expected to start when C = Vmax = 1- Using 
this value for C we obtain from I 



^yy 



therefore 6 = 6 fVA'^+i g^j^^ from [B2| we then have 



(B3) 



K - r[-VK^ + 1 + 1] - G(l - r) = 



(B4) 



which by solving gives us 



G(l - r) + r - r^{G{l-r)+ry + (1 - r^) 



l-r2 



Numerical integration confirms this result. 



(B5) 
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APPENDIX C: KH-MODES IN THE SMALL pi/p2 LIMIT 



We begin with the Rayleigh equation 3.3 for large K 



'yy 



Set z = Ky and e = 1 //C ; we then have 

V.Jez) 



1 + e 



2 ^ zz\ 



V{ez)-C^ 



or 



<p..- {l+e^F{ez)) = 0, 
where F{x) = Vzz{x)/{V{x) — C). Expanding (f) in powers of 



i)o + eVi + £^2 + - , 



we obtain, to first order, 



^-0 = e-^ ; 



at the next order, we have 



'Izz — VI 



F{ez) 



which has the solution 

01 = y G{\z - x\)F{tx)(t)Q{x)dx + Ae'"" , 

where G = — ^e"'"^"^' is the Greens function and A is a constant to be chosen so that (pi will 
satisfy the boundary condition at zero: 



1 r e-(^-")e-"F(ex)rfx - l T e~^'-'-'^e~^F{ex)dx + Ae 

2 Jo 2 Jz 



\ rz 1 roo /I roo 

--e-' / F{ex)dx - -e' / e-^''F{ex)dx + - / e-^''F{ex)dx \ e 

2 Jo 2 Jz \2 Jo 
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— e 



rz poo poo 

/ F{ex)dx+ / e-'^'^F{e{w + z))dw - / e''^'' F{ex)dx 
Jo Jo Jo 



{w = X — z) 



We are interested in the first derivative of at zero, so we can obtain 



dz 



z=0 



1 _^d_ 
'2^ d^ 



pz poo poo 

/ F{ex)dx + / e-^'"F{e{w + z))dw - / e''^'' F{ex)dx 
Jo Jo Jo 



2 = 



2' 



F{0) + j-(e^' e-^^F{ex)dx 
dz ^ J z 



2=0 



F(0) + 2 / e-'^''F{ex)dx - F(0) 
Jo 



2=0 



roc 

^0 



2 = 



/f e-^(F(0) + e|F'(0) +e^l {^J F"iO) + ...)d 



'2 



.1f(0)- V(0) + ie^F"(0) + ... 



Tlierefore tlie final result for the first derivative of (j) at zero is 



^2/ 1 2=0 



2V0-C 4 



oy 10 



+ ... 



Applying the boundary condition 3.4 at the interface 



KC^ - r[{Vo - Cf(f)y\o - (Vo - c)V'\o] - G(l - r) = 



we obtain 



KC^ _ r[-K{Vo - Cy - K-'^{K - C)V"\o + ... - - C)V'\o] - G(l - r) = 
Scaling K and C so that K = k/r and C = v^c, and substituting we have 
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kc'' + k{Vo - V^cf + r^k-^hVo - V^c)V"\o + r{Vo - V^c)V'\o - G + rG = 
If we expand c in powers of r^/^, 

c = Co + ci/2r^/^ + cir + Cg/sr^/^ + ... , 
we can obtain each term separately. Here we give only the first few terms: 



Co: k4 + kV„'^G^Co^ JG/k-V^^ ; 



ci/2 : 2kcoCi/2 - 2kVoCo = ^ ci/2 = K ; 



ci : A;c^/2 + 2A;coCi + kcl - 2A;Kci/2 + kc^ + VoV'\o + G = 



3 1 

ci = --Co - -v;y'|o/co 
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(a) 



G = 0.2 
G = 
G = 
= 
0=1 










FIG. 4. Imaginary part of C for a logarithmic wind profile and V^=1.0 (a)r = 0.01, (b)r = 0.1, 
(c)r = 0.5 . 
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FIG. 5. Imaginary part of C for a logarithmic wind profile with 

(b)yo=o.5, (c)yo=i.o, (d) ^0=1.5 . 
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FIG. 7. Growth rate of the fastest growing mode as a function of G (a) logarithmic wind profile, 
(b) tanh wind profile . 




FIG. 8. Slopes of the previous graph as a function of r (e.g. the dependence of the exponent 
in eq. 3.15| on r). 
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FIG. 9. Intersections of the previous graph as a function of r (e.g. the dependence of the 
coefficient in front the exponential in eq. |3.15| , on r). 
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FIG. 10. The dependence of the wavenumber of the fastest growing mode on G (a) logarithmic 
wind profile, (b) tanh wind profile. 
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FIG. 11. Imaginary part of C for a logarithmic wind profile with G=0.5 in the presence 
surface tension (a)r = 0.01, (b)r = 0.1, (c)r = 0.5 . 
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FIG. 12. Imaginary part of C for a tanh wind profile with G=0.5 in the presence of surface 
tension (a)r = 0.01, (b)r = 0.1, (c)r = 0.5 . 
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FIG. 13. Im{C} for r=0.1 for the compressible case (a)G=0.1, (b)G=0.01 
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TABLES 



TABLE I. Approximate range for parameter G in three different situations. 

Ui(cm s~^) g cm s~^ a~^cm G 

ocean 10^ - 10^ 10^ 10 - 10^ IQ-^ ~ 1 

Sun's surface 10^ ~ 10^ 10^ =^ 10^ ~ 10^ 10^ =^ ~ lO^-^ 

WD 10^ ~ 10^ 10^ 10=* ~ 10^ 1 ~ 10 



45 



